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ABSTRACT 

The  aims  of  this  paper  are  threefold:  to  increase  the  level  of  awareness  within  the 
shock  capturing  community  to  the  fact  that  many  Godunov-type  methods  contain 
subtle  flaws  that  can  cause  spurious  solutions  to  be  computed;  to  identify  one  mech¬ 
anism  that  might  thwart  attempts  to  produce  very  high  resolution  simulations;  and 
to  proffer  a  simple  strategy  for  overcoming  the  specific  failings  of  individual  Riemann 
solvers. 
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1  Introduction 

Over  recent  years,  a  plethora  of  shock-capturing  schemes  have  been  developed  for 
the  Euler  equations  of  gas  dynamics.  During  this  period,  it  has  emerged  that  one 
of  the  more  successful  strategies  for  designing  a  shock  capturing  scheme  is  to  follow 
Godunov’s[19]  lead  and  utilize  a  classic  initial-value  problem  known  as  a  Riemann 
problem[3].  Godunov  assumed  that  a  flow  solution  could  be  represented  by  a  series  of 
piecewise  constant  states.  Thus,  the  numerical  representation  closely  approximates 
the  true  solution  near  discontinuities,  and  regions  of  smooth  flow  are  reasonably  well 
approximated  by  a  series  of  step  functions.  He  evolved  this  discretized  flow  solution 
by  considering  the  nonlinear  interactions  between  its  component  states.  Viewed  in 
isolation,  each  pair  of  neighbouring  states  constitutes  a  Riemann  problem,  the  .solution 
to  which  may  be  found  exactly [8].  The  results  from  these  separate  Riemann  problems 
may  then  be  averaged  so  as  to  advance  the  flow  solution  through  ome  time  increment. 
Because  it  mimics  much  of  the  relevant  physics,  Godunov’s  scheme  results  in  an 
accurate  and  well-behaved  treatment  of  shock  waves. 

Although  it  provides  the  bedrock  upon  which  most  modern  schemes  are  built,  in 
its  original  form  Godunov’s  method  is  o'’  limited  use.  Firstly,  the  scheme  proves  to 
be  highly  dissipative  and  so  it  requires  an  inordinately  fine  mesh  to  resolve  complex 
shock-on-shock  interactions.  Secondly,  since  a  Riemann  problem  has  no  closed  form 
solution  and  can  only  be  solved  by  some  iterative  method,  Godunov’s  scheme  is 
significantly  more  expensive  than  schemes  which  employ  ordinary  finite-difference 
operators. 

One  of  the  first  people  to  address  this  second  shortcoming  was  Roe[19].  He  argued 
that  since  the  Riemann  problems  associated  with  Godunov’s  method  arise  from  an 
approximation  of  the  data,  it  might  be  sufficient  to  find  only  approximate  solutions 
to  these  Riemann  problems,  provided  that  they  still  describe  important,  nonlinear 
behaviour.  His  motivation  being,  that  approximate  solutions  can  be  computed  much 
more  cheaply  than  exact  solutions.  Thus  the  industry  of  designing  approximate  Rie¬ 
mann  solvers  was  born[14,  15,  4,  23].  Now  whilst  Godunov-type  schemes  are  often 
held  up  to  be  models  of  robustness  they  can  on  occasions  fail  quite  spectacularly.  For 
example,  when  computing  shock  reflection  problems.  Roe’s  method  can  sometimes  go 
awry  by  admitting  solutions  for  which  the  Mach  stem  is  inexplicably  kinked.  The  ex¬ 
istence  of  such  failings  partly  explains  why  no  consensus  of  opinion  has  been  reached 
concerning  the  ideal  Ri<"mann  =olvcr.  Whenever  a  new  failing  is  unearthed,  it  adds 
fuel  to  the  great  Riemann  solver  debate;  method  ,V  is  better  than  method  V’.  because 


of  reasons  A,  B  and  C.  It  is  our  contention  that  for  all  the  current  crop  of  Rieinaiin 
solvers,  at  least  one  set  of  circumstances  may  be  fonnd  for  which  the  solver  is  found 
wanting;  some  failings  are  just  more  visible  than  others. 

In  Section  2  we  catalogue  a  number  of  situations  in  which  anomalous  behaviour  is 
known  to  occur.  This  catalogue  should  serve  to  increase  the  general  levfd  of  awareiu'ss 
within  the  shock  capturing  community  to  the  current  limitatioTis  of  Riemann  solver 
technology.  At  present,  this  awareness  is  not  all  that  it  shoidd  be,  there  are  many 
instances  in  the  literature  where  suspect  numerical  results  are  presented  with  either 
little  or  no  adverse  comment.  VVe  believe  that  one  of  the  failings  listed  m  our  catalogue 
has  hitherto  gone  unreported.  In  Section  2  we  proffer  a  diagnosis  of  the  mechanism 
wtiich  caii'-'Ms  this  new  failing. 

At  this  juncture  it  should  be  noted  that  any  foibles  that  a  specific  Riemann  solver 
might  have,  may  usually  be  controlled  by  the  judicious  use  of  a  small  amount  of 
artificial  dissipation.  Indeed,  it  is  worth  pointing  out  that  Woodward  and  Collela’s 
PPM  scheme[2],  which  has  proved  itself  to  be  more  robust  than  most  Codunov-type 
schemes,  does  in  fact  employ  an  elaborate  artificial  dissipation  mechanism  to  supple¬ 
ment  the  dissipation  provided  via  upwinding.  As  will  be  described  in  Section  4,  we 
favour  a  strategy  whereby  the  weaknesses  of  any  one  solver  are  overcome  by  combin¬ 
ing  it  with  one  or  more  complementary  solvers.  The  main  advantages  of  this  approach 
over  that  of  adding  artificial  dissipation  arc  twofold.  Firstly,  it  does  not  degrade  the 
resolution  of  the  base  Riemann  solver;  it  is  possible  to  control  certain  instabilities  by 
changing  the  flavour  of  the  dissipation  mechanism  rather  than  increasing  the  absolute 
level  of  dissipation.  Secondly,  it  does  not  necessitate  a  host  of  tunable  parameters, 
and  so  this  synergetic  strategy  does  not  negate  the  principal  advantage  of  Godunov- 
type  schemes  over  other  shock  capturing  methods.  Of  course,  we  are  left  with  the 
difficulty  of  deciding  when  to  use  one  Riemann  solver  in  preference  to  another;  how¬ 
ever,  we  present  a  number  of  computations  which  suggest  that  this  difficulty  is  not 
particularly  bothersome. 

Finally,  in  Section  5,  we  list  the  main  conclusicns  that  we  have  drawn  from  this 
work.  Note  that  in  this  paper  we  do  not  address  the  first  sliortcoming  of  Godunov's 
method,  namely  its  low  resolution.  Following  van  Leer[26],  it  is  assumed  that  a  high- 
order  extension  to  a  first-order  method  can  always  be  achieved  by  pre-processing  the 
data  supplied  to  the  Riemann  .solver. 


2  A  Catalogue  of  Failings 
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Wc  now  present  several  instances  where  various  Rieniann  solvers  ar(’  known  to  giv(' 
unreliable  results.  While  most  of  these  cases  will  be  known  to  the  aficionado,  we 
believe  that  one  of  the  cases  which  we  are  about  to  describe  has  hitherto  gone  un¬ 
reported  in  the  literature.  Although  our  catalogue  is  not  exhaustive,  we  hope  that 
it  might  save  some  investigators  from  the  harrowing  experience  of  spending  weeks  or 
even  months  searching  for  coding  errors  that  simply  do  not  exist. 

All  the  computational  results  shown  in  this  section  are  for  first-order  schemes, 
since  such  methods  have  low  resolution  our  calculations  employed  relatively  fine 
meshes;  for  clarity,  grids  are  drawn  using  every  other  fourth  grid  line. 

2.1  Expansion  Shocks 

By  far  the  most  widely  investigated  failing  is  that  some  Riemann  solvers  do  not 
satisfy  an  “entropy  condition”  such  schemes  can  admit  non-physical  solutions  such 
as  expansion  shocks.  0sher[14]  has  found  a  general  condition  for  a  scheme  to  be 
entropy  satisfying  when  applied  to  scalar  equations  and  he  designates  such  schemes  E- 
schemes.  At  present,  however,  any  extension  to  a  system  of  equations  contains  a  large 
amount  of  empiricism  and  must  therefore  remain  suspect.  Indeed,  Godunov’s  method 
is  classified  as  an  E-scheme  but,  as  observed  by  Woodward  and  Collela[27],  it  can  give 
rise  to  nearly  discontinuous  expansion  fans  near  sonic  points.  The  density  contours 
shown  in  Figure  1  illustrate  this  deficiency  of  Godunov’s  method  quite  clearly.  These 
results  are  taken  from  the  diffraction  of  a  strong  shock  wave,  Mg  —  5.09  with  7  =  1.4, 
around  a  90°  corner. 

In  its  basic  form.  Roe’s  scheme  is  another  solver  that  admits  expansion  shocks, 
however,  several  fixes  have  been  proffered  which  cure  Roe’s  scheme  of  this  particu¬ 
lar  affliction [20,  5,  28].  Such  fixes  are  typical  of  the  way  in  which  Riemann-solver 
deficiencies  have  been  countered  up  to  now.  Whilst  this  strategy  has  proved  reason¬ 
ably  successful,  it  has  a  number  of  drawbacks.  Sometimes  a  fix  uses  a  parameter 
which  must  be  retuned  between  problems  and  hence  one  of  the  major  advantages  of 
Riemann-based  schemes  over  say  artificial-dissipation  schemes  is  lost.  Alternatively, 
a  scheme  may  require  more  than  one  fix,  and  it  may  be  unclear  how  the  different  fixes 
interact  with  one  another. 


(a)  (b) 


Figure  1:  A  strong  shock  diffracting  around  a  corner  gives  rise  to  an  expansion  shock: 
(a)  Density  contours,  (b)  Computational  grid. 

2.2  Negative  Internal  Energies 

Another  situation,  in  which  some  Riemann  solvers  are  found  wanting,  occurs  when¬ 
ever  the  dominant  energy  mode  is  kinetic  rather  than  thermal.  For  such  solvers, 
the  kinetic  energy  computed  from  a  numerical  approximation  to  the  conservation 
laws  of  mass  and  momentum,  can  exceed  the  total  energy  computed  via  an  approx¬ 
imation  to  the  conservation  law  of  energy.  Thus  they  can  yield  negative  internal 
energies,  and  hence  negative  pressures,  which  cause  the  scheme  to  fail.  Einfeldt  ct 
a/.[5]  call  any  scheme  which  can  be  guaranteed  not  to  yield  negative  pressures,  “pos¬ 
itively  conservative”.  They  have  shown  that  while  Godunov’s  scheme  is  “positively 
conservative”,  the  reverse  is  true  for  any  Godunov-type  scheme  based  on  a  linearized 
Riemann  solver.  Indeed,  the  basic  form  of  Roe’s  scheme  is  unable  to  cope  with  the 
test  problem  shown  in  Figure  1;  the  strength  of  the  diffracting  shock  is  sufficient  to 
cause  a  negative  pressure  to  be  computed  near  the  apex  of  the  corner.  Roe’s  scheme 
may  be  made  “positively  conservative” by  modifying  its  wave  speeds,  in  essence,  the 
scheme  is  made  more  dissipative  by  increasing  the  spread  in  velocity  between  the  two 
acoustic  waves[5]. 


2.3  Slowly  Moving  Shocks 

Since  shock  capturing  schemes  do  not  resolve  the  internal  structure  of  a  shock  wave, 
no  physical  significance  can  be  attached  to  the  discrete  shock  structure  produced  by  a 
numerical  scheme;  methods  are  built  upon  the  premise  that  shock  profiles  are  mono¬ 
tone,  the  precise  structure  comes  about  as  a  matter  of  course  and  is  not  preordained. 
Unfortunately,  Roberts  has  shown  that  the  nature  of  the  shock  structure  produced 
by  a  particular  scheme  can  have  a  large  bearing  on  how  well  the  scheme  copes  with 
slowly  moving  shock  waves[21].  Godunov-type  methods  fare  quite  badly  in  this  re¬ 
spect,  as  the  shock  moves  relative  to  the  mesh,  the  shock  profile  flexes,  perturbing 
the  supposedly  passive  characteristic  fields  as  it  does  so. 

Figure  2  shows  a  snapshot  of  the  shock  profile  produced  by  Einfeldt’s  HLLE 
sch^me[4],  taken  from  the  simulation  of  a  shock  wave  which  is  moving  slowly  from 
left  to  right;  the  pre-shock  state,  (density,  velocity,  pressure),  is  (1,— 3.44,1),  and 
the  post-shock  state  is  (3.86,-0.81,  10.33).  Note  that  for  a  Courant  number  of  one, 
it  takes  50  time  steps  for  this  shock  to  traverse  one  mesh  cell.  The  low  frequency 
perturbations  observed  in  this  figure  are  produced  to  a  greater  or  lesser  extent  by 
any  scheme  which  attempts  to  “iecognize”a  shock  wave.  For  fast  moving  shocks, 
the  post-shock  noise  will  be  of  a  much  shorter  wavelength  than  is  the  case  here,  and 
will  be  effectively  damped  by  the  dissipation  of  the  scheme.  Roberts  reports  that 
Osher’s  <:chcme[I3]  does  not  produce  low  frequency  noise  for  slowly  moving  shocks, 
since  it  never  connects  two  adjacent  states  by  a  shock,  and  he  concludes  that  there 
may  be  advantages  to  using  flux  formulas  that  do  not  recognize  the  analytic  shock 
jump  conditions. 

Another  situation,  where  the  perturbation  of  a  shock  from  its  preferred  profile 
results  in  perturbations  on  the  passive  characteristic  fields,  occurs  whenever  a  shock 
crosses  a  discontinuity  in  mesh  spacing[17].  But,  in  this  case,  sizeable  perturbations 
may  occur  whatever  the  speed  of  the  shock. 


Figure  2;  Low  frequency,  post-shock  oscillations  occur  for  slowly  moving  shock  waves. 
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2.4  The  Carbuncle  Phenomenon 

Several  authors  have  by  now  reported  a  failing  of  Roe’s  scheme  which  has  been  dubbed 
the  “carbuncle  phenomenon” [16,  12,  11].  For  steady-state,  blunt-body  calculations, 
Roe’s  scheme  sometimes  admits  a  spurious  solution  in  which  a  protuberance  grows 
ahead  of  the  bow  shock,  along  the  stagnation  line.  It  appears  that  this  effect  is 
more  pronounced  the  more  closely  the  grid  is  aligned  to  the  bow  shock.  Also,  a 
carbuncle  is  more  likely  to  appear  for  high  Mach  number  flows  than  for  low  Mach 
number  flows.  Figure  3  shows  such  a  spurious  solution,  here  the  freestream  Mach 
number  was  taken  to  be  10.  Note  that  along  the  stagnation  line,  the  bow  shock  is 
almost  perfectly  aligned  with  the  grid.  Consequently,  parallel  to  the  shock.  Roe’s 
scheme  will  not  add  any  dissipation  via  the  contact  and  shear  waves,  to  counteract 
perturbations  that  appea"  through  the  acoustic  waves;  this  appears  to  be  a  recurring 
theme  whenever  Roe’s  method  fails.  It  is  interesting  to  note  that  if  Marten’s  entropy 
fix[28]  is  applied  to  the  contact  and  shear  waves,  any  shortcoming  of  Roe’s  scheme  is 
invariably  cured.  However,  there  is  no  justification,  either  physical  or  mathematical, 
for  applying  this  fix  to  these  waves,  it  is  just  a  convenient  method  for  introducing  an 
amount  of  artificial  dissipation  into  the  scheme. 

2.5  Kinked  Mach  Stems 

During  the  course  of  developing  a  mesh  adaption  scheme,  we  encountered  a  failing  of 
Roe’s  scheme  which  is  not  dissimilar  to  the  “carbuncle  phenomenon”[17].  When  the 
reflection  of  a  plane  shock  wave  from  a  ramp  lies  in  the  double  Mach  reflection  regime, 
the  principal  Mach  stem  is  sometimes  inexplicably  kinked.  Figure  4  shows  a  snapshot 
of  the  pressure  contours  taken  during  the  reflection  of  a  plane  shock,  Ms  =  5.5  with 
7  =  1.4,  from  a  30"  ramp.  The  principal  Mach  stem  is  so  severely  kinked  that  it  has 
given  rise  to  a  spurious  triple  point.  Similarly  strange  results  have  been  produced  by 
Sawada[22],  and  by  Itoh  and  Takayama[9].  As  before,  because  of  the  way  the  Mach 
stem  is  aligned  with  the  grid,  there  is  probably  insufficient  dissipation  added  via  the 
contact  and  shear  waves  to  counteract  perturbations  that  appear  via  the  acoustic 


waves. 


Figure  4:  The  principal  Mach  stem  arising  from  the  refli'ction  of  plane  shock  from 
ramp  is  inexplicably  kinked:  (a)  Pressure  contours,  (b)  Computational  grid. 
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2.6  Odd-Even  Decoupling 

By  far  the  most  insidious  failing  that  we  have  come  across  has.  we  helicvc.  none 
unreported  in  the  literature.  During  the  course  of  piaxlucing  very  liigli  resolution 
simulations,  we  have  noticed  a  teiuh'ncv  for  odd-('v«’n  decoupling  to  occut  along  the 
length  of  planar  shocks  which  are  aligiuxl  with  the  mesh,  for  an  example  six*  l-'igure  1  1 . 
Of  the  Riemann  solvers  that  we  ha\e  at  our  tlispo.sal.  this  failing  alllicts;  an  ('xact 
solver[21].  Roe's  solver  and  Toro's  rmeariz<'d  solver[‘2d].  VVe  ('lupliasizf'  the  fact  that 
this  phenomena  only  hecumes  apparent  for  very  high  resolution  simulations  which 
suffer  .some  systematic  pertuhation.  However,  as  will  he  shown  Ix'low.  tin-  requin'd 
perturbation  can  arise  quite  innocuously,  .so  we  suspect  that  this  failing  will  prov*' 
fairly  widespread  once  very  high  resolution  simulatioiis  heconu'  comtnonplace  due  to 
increa.ses  in  computer  power. 

.Now  since  we  obtain  our  high  grid  resolution  by  mt'ans  of  a  fairly  comph'X  mesh 
adaption  .scheme[17].  it  seemed  rea.sonable  to  suppo.se  that  this  odd-cnen  decoupling 
was  attributable  to  some  coding  error,  but  an  exhaustive  search  for  such  an  error 
proved  fruitless.  Subsec|uently,  we  have  managed  to  reproduce  this  failitig  in  a  more 
controlled  manner,  as  has  a  colleague  using  an  independent  code[in].  so  we  have  little 
doubt  that  this  tendency  for  odd-even  decoupling  to  occur,  constitutes  a  genuine 
failing,  rather  than  being  the  manifestation  of  some  deficiency  of  our  code,  fhat 
said,  our  arlaptive  mesh  scheme  clearly  exasperate;;  this  failing.  In  the  tu'xt  si'ction 
we  shall  present  a  possible  diagnosis  of  the  mechanism  which  cause's  this  mode'  of 
failure,  here  we  merely  present  the  evidence  that  it  exists. 

Figure  5  shows  several  snapshots  of  the  density  contours  from  the  simulation  of 
a  plane  shock  wave,  .V/j  =  6  with  y  =  1.1.  propagating  down  a  fluct.  For  this 
calculation  v.?  have  used  Roe’s  scheme,  together  with  a  nominally  uniform  grid  of  20 
by  800  cells,  with  unit  spacing,  the  centre-line  of  which  is  perturbed  from  that  of  a 
perfectly  uniform  mesh  in  the  following  manner 

^ jmid  T  10  *  for  /  eVl'Il. 

=  y^mid  —  10  for  i  odd. 

This  perturbation  to  the  grid  centre-line  promotes  odd-even  decoupling  along  the 
length  of  the  shock.  Note  that  the  shock  has  projiagated  some  ITi  channel  widths 
before  the  decoupling  first  becomes  apparent,  see  frame  (b).  For  this  point  in  the 
calodation,  Figure  6  shows  a  series  of  slices  across  the  duct,  for  both  the  dc'usity  and 
pressure  fields,  as  one  moves  from  the  head  to  the  foot  of  the  shock.  Intt'rest  ingly. 


within  the  shuck  tlie  decuiiplinfj;  of  tlu'  pressure  field  is  out  of  phase  with  the  decou 
pliiig  of  the  dettsity  field.  .\s  the  shock  continues  to  |»ropagatc  (h)W!;  the  duct.  s( 
the  dtMoupling  Ix'coiues  propo-.dvt'ly  worse  until  the-  shock  hia'aks  down  coinj)letely 
However,  at  no  stage  in  *'  .  calculation  does  the  codi'  blow  uj)  in  the  st'iisc  that  it 
geiK'rates  a  floating  point  exc('ption:  it  siin|)Iv  go<’s  astray. 


Odd-even  dt'coupling  occurs  for  a  shock  propagating  down  a  dm  t 


Cell  Across  Duct 


3  Odd-Even  Decoupling  -  A  Diagnosis? 


As  ye't,  tlu'  tools  ilo  not  exist  that  would  allow  us  to  perforin  a  rigorous  nonlinear 
stability  analysis  for  some  shock  capturing  scheme  applied  to  the  Euler  et|uations. 
However,  it  is  possible  to  examine  the  way  in  which  a  scheme  evolves  certain  sets 
of  prescribed  data,  so  as  to  ascertain  its  likely  stability  characteristics.  Since  some 
Riemann  solvers  allow  odd-even  decoupling  to  develop  along  the  length  of  a  plane 
shock,  it  might  prove  fruitful  to  examine  how  different  schemes  evolve  sawtooth-type 
initial  tlata.  To  this  end,  we  consider  the  one-dimensional  Euler  ecpiations  with  a 
passive  component  of  shear  velocity 


dt 


(  pv  \ 

pa 

d 

^dy 

puv 

pv 

px'^  +  P 

V  E  ] 

V  (£■  +  P)v  ! 

=  0. 


(1) 


The  cptantities  p.  p.  u.  v  and  E  are  density,  pressure,  the  passive  shear  component 
of  velocity,  the  velocity  in  the  y  direction,  and  the  total  energy  per  unit  volume, 
respectively.  For  a  perfect  gas 

E  = 

—  1  ^ 


where  -y  is  the  ratio  of  specific  heats.  We  assume  that  the  computational  mesh  is 
uniform,  w'ith  mesh  spacing  iA//,  and  that  the  discrete  solution  at  time  is  given  by 
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if  j  is  even,  and 
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if  j  is  odd.  Here  /5"  and  are  the  amplitudes  of  the  sawtooth  profiles  for  the  density 
and  the  pressure  fields.  We  shall  consider  two  schemes  which  may  be  expressed  in 
the  form 


w;"  ’  =  w" 


St 

■A?/ 


(3) 


where  W  is  the  con.served  variable  vector  {p,pxi,px\  Ey.  and  G’‘  ^  is  a  first-order  flux 
function  romj)Ute(l  from  tlie  states  W"  and  W"^,. 
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3.1  Roe’s  scheme 


Using  Roe’s  scheme[19],  the  interface  flux  for  the  system  of  equations  (1)  may  he 
written 

fc=4 


k=\  ^ 


where  the  wave  speeds  {AU)},  wave  strengths  {oUl},  and  eigenvectors  {eU>}  are  given 

by 

Ab)  z=  V  —  a,  A*^^  =  AUI  =  u,  A*'^'  =  e  +  d; 
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Here,  quantities  written  as  (•)  are  the  so-called  Roe  averaged  quantities,  a  and  k  are 
the  Roe  averaged  sound  speed  and  total  enthalpy  respectively,  and  A(»)  represents 
the  forward  difference  operator  —  (•)j.  Now  for  our  chosen  data. 


and  A(*).^i  = 


Therefore, 


k=4 

E 

fc=i 


yW 

J  +  5 


k=4 

eW  _  Jk) 


J  2 


e'tt. 

2 


Also,  G"_,  =  G”^i,  so  the  evolution  scheme  (3)  may  be  written 
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which  can  simplified  to 
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Where,  Uy  is  the  Courant  number  — — .  Recognizing  that  W"  may  be  expressed  as 
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and  that  by  definition 
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equation  (4)  may  be  manipulated  to  give, 


'n+l  *»i 
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and 

r*'  =  -  2,/,). 

From  which  it  can  be  seen  that  the  initial  perturbation  to  the  pressure  field  is  damped, 
provided  that  the  CFL  condition  is  met,  that  is 


Vy  <  1. 

However,  the  form  of  the  evolution  for  the  density  perturbation  exposes  a  flaw  in 
Roe’s  scheme;  the  density  perturbation  is  fed  directly  from  the  pressure  perturbation. 
Making  the  loose  approximation  that  a  remains  constant,  for  a  one-off  disturbance 
{p\p^) 


- j^l  -f  ( 1  —  2vy)  +  ( 1  —  -I-  ( 1 


2;/'y) 


Tl-l 


Thus, 


p- = p"  -  C- 

Therefore,  if  p°  is  of  opposite  sign  to  then  for  a  one-off  disturbance  p  will  grow 
but  it  will  remain  bounded.  But  if  the  pressure  field  is  continuously  perturbed  in 
a  systematic  manner,  no  matter  how  small  the  pressure  perturbations,  p  will  grow 
without  bound,  albeit  slowly.  For  two-dimensional  calculations,  although  we  cannot 


prove  it,  we  suspect  that  a  strong  shock  wave  moving  normal  to  the  y  direction 
provides  this  systematic  perturbation. 

Firstly,  it  is  interesting  to  note  that  the  failing  reported  in  Section  2.6  is  only 
observed  for  strong  shocks.  For  a  strong  shock  wave,  it  seems  reasonable  to  suppose 
that  II2 1  is  more  likely  to  be  larger  than  |p®|,  than  would  be  the  case  for  weak  shocks. 
So,  even  if  and  fp  are  initially  of  the  same  sign,  they  need  not  remain  so,  see  (5). 
Now  consider  frame  (b)  of  Figure  5  and  the  associated  profiles  shown  in  Figure  6, 
within  the  shock  the  odd-even  decoupling  of  the  pressure  and  density  fields  are  indeed 
out  of  phase  with  one  another,  which  is  consistent  with  the  observations  made  above. 
Such  behaviour  will  cause  the  local  sound  speed  to  vary  along  the  length  of  the 
shock,  its  profile  will  exhibit  a  sawtooth  perturbation  which  is  in  phase  with  that  of 
the  pressure  field.  Consequently,  the  individual  segments  of  the  shock  will  be  moving 
alternately  faster  and  slower  than  the  nominal  shock  speed.  Such  movements  will 
exaggerate  the  sawtooth  perturbation  to  the  pressure  field  along  the  length  of  the 
shock,  but  diminish  that  for  the  density  field.  The  increased  pressure  perturbations 
will  then  promote  an  increase  in  the  density  perturbations  as  detailed  above,  and  so 
the  whole  process  repeats  itself. 

Since  there  are  two  competing  processes  that  affect  the  density  perturbations, 
namely,  the  relative  movements  of  the  shock  and  the  decoupling  along  the  length  of 
the  shock,  we  cannot  categorically  state  that  Roe’s  method  is  bound  to  break  down. 
However,  the  weight  of  numerical  evidence  suggests,  that  at  least  for  strong  shocks 
Roe’s  scheme  will  break  down  in  the  manner  described  here.  Given  our  arguments,  it 
should  come  as  no  surprise  that  Godunov’s  method  also  exhibits  a  tendency  to  allow 
odd-even  decoupling  to  occur  along  the  length  of  a  stron^,  shock  wave.  Since  it  is 
the  sweep  parallel  to  the  shock  that  primarily  causes  the  instability,  the  differences 
between  using  an  exact  Riemann  solver  and  Roe’s  linearized  solver  for  data  that  is 
nominally  uniform  should  have  little  bearing  on  the  growth  of  the  instability. 

Finally,  before  moving  on,  it  should  be  noted  that  none  of  the  popular  entropy 
fixes  which  are  applied  to  Roe’s  scheme  cure  this  particular  failing,  excepting  the 
case  where  Marten’s  fix  is  applied  to  the  shear  and  contact  waves^;  simply  altering 
the  acoustic  wave  speeds  can  have  no  affect,  because  of  the  symmetry  of  the  data,  both 
waves  will  be  changed  by  the  same  amount,  and  so  the  problem  will  persist.  Also, 
moving  to  a  high-order  version  of  Roe’s  scheme  will  not  improve  matters,  because  the 
odd-even  decoupling  will  cause  a  high-order  flux  function  to  drop  to  the  first-order 
function. 

^To  reiterate  the  comment  made  in  .Section  2.4,  applying  Marten’s  entropy  fix  to  the  linearly 
degenerate  wave  fields  has  no  mathematical  or  physical  justification,  it  is  merely  a  convenient  way- 
in  which  to  add  an  amount  of  artificial  dissipation. 
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3.2  Einfeldt’s  HLLE  scheme 

For  Einfeldt’s  HLLE  scheme[4],  the  interface  flux  function  is  given  by 


G"  1  = - — - - - fW"  -W’M 

6(+)-6{-) 


where 
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and  A('),  At**)  are  the  two  acoustic  wave  speeds  from  Roe’s  method.  Now  for  our 
chosen  data,  (•), .  j.  is  equal  to  (•)  _!,  therefore 
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Using  these  signal  weightings,  it  may  be  found  that 
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From  which  it  can  be  seen  that  both  the  density  and  pressure  perturbations  are 
damped,  provided  that  the  CFL  condition  is  satisfied.  Just  as  important,  however,  is 
the  fact  that  the  pressure  perturbation  does  not  feed  into  the  density  perturbation,  so 
we  would  not  expect  the  HLLE  solver  to  exhibit  the  odd-even  decoupling  that  afflicts 
both  Roe’s  scheme  and  Godunov’s  scheme;  numerical  experimentation  confirms  this 
expectation. 

It  is  our  contention  that  any  scheme  for  which  it  can  be  shown  that  the  perturba¬ 
tion  to  the  pressure  field  feeds  the  perturbation  to  the  density  field,  will  be  afflicted 
by  the  odd-even  decoupling  shown  in  Figure  5.  Thus  it  comes  as  no  surprise  to  find 
that  Toro’s  linearized  Riemann  solver[23]  is  afflicted  by  this  failing,  but  Lion’s  AUSM 
scheme[ll]  is  not.  The  way  is  now  open  for  some  mathematician  to  perform  a  more 
rigorous  analysis  than  we  are  able,  so  as  to  shed  additional  light  on  the  mechanism 
which  causes  this  particular  failing. 
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4  An  Adaptive  Riemann  Solver 

Having  exposed  many  of  the  weaknesses  of  Riemann  solvers,  we  now  present  a  simple 
strategy,  that  we  have  found  useful,  for  improving  the  robustness  of  Godunov-type 
schemes.  In  essence,  we  select  the  precise  flavour  of  upwinding  to  match  the  local 
flow  data,  such  that  a  particular  Riemann  solver  is  only  employed  in  those  situations 
where  it  is  known  to  give  reliable  results.  By  recognizing  the  limitations  of  any  one 
solver  it  is  possible  to  reap  its  advantages  without  suffering  its  attendant  failings. 

Our  synergetic  strategy  has  a  number  of  attractions,  not  least  of  which  is  that  some 
favoured  solver  need  not  be  jettisoned  simply  because  it,  occasionally,  fails.  However, 
it  does  introduce  the  difficulty  of  how  to  decide  when  to  use  one  Riemann  solver  in 
preference  to  another.  But  it  has  been  our  experience  that  this  added  difficulty  is  not 
particularly  bothersome,  for  we  tend  to  combine  a  single  high  resolution  Riemann 
solver  with  just  one  or  two  other  solvers  that  prove  more  reliable  under  conditions 
which  are  fairly  well-defined,  and  so  a  set  of  ad  hoc  switching  functions  suffice.  For 
example,  some  of  the  worst  failings  of  Riemann  solvers  occur  in  the  vicinity  of  strong 
shock  waves.  To  overcome  such  failings  we  employ  Einfeldt’s  HLLE  scheme.  Now  it 
makes  little  sense  to  chop  and  change  the  choice  of  Riemann  solver  used  along  the 
length  of  a  shock  wave,  since  to  do  so  would  inevitably  perturb  a  planar  shock  front. 
Hence,  we  apply  this  particular  Riemann  solver  throughout  the  immediate  vicinity  of 
a  strong  shock.  Thus  the  HLLE  switching  function  need  only  locate  the  position  of  a 
shock  wave,  but  such  functions  already  exist  in  the  guise  of  mesh  refinement,  monitor 
functions. 

A  simple  test  that  identifies  those  cell  interfaces  which  are  in  the  vicinity  of  a 
strong  shock  is  to  check  whether  or  not 


\Pr  -  Pi\ 
min{pi,pr) 


> 


(6) 


where  a  is  some  threshold  parameter  which  is  problem  dependent  and  pr  and  pi  refer 
to  the  pressures  which  act  on  the  interface.  If  this  condition  is  met,  the  two  cells 
separated  by  the  interface  are  flagged  a.s  lying  within  a  strong  shock.  So,  when  it 
comes  to  computing  cell-interface  fluxes,  if  the  cells  either  side  of  an  interface  are 
both  flagged  as  lying  within  a  strong  shock,  the  flux  is  computed  using  the  HLLE 
solver.  Note  that  since  numerical  shocks  are  invariably  smeared  over  several  mesh 
cells,  it  is  worth  locating  shocks  using  a  projection  of  the  flow  solution  on  a  grid  which 
is  coarser  than  that  used  for  the  calculation.  On  such  a  grid  a  shock  will  appear  much 
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less  smeared,  and  so  the  left-hand  side  of  the  above  switching  function  will  be  a  fair 
indication  of  its  strength.  Once  a  set  of  cells  have  been  flagged  on  this  coarse  mesh, 
the  flags  may  be  prolongated  to  the  actual  computational  mesh  so  as  to  find  those 
cells  which  lie  in  the  vicinity  of  a  shock  wave. 

Before  proceeding  further,  several  observations  should  be  made.  Firstly,  although 
the  HLLE  solver  is  adjudged  to  be  a  low  resolution  scheme,  it  does  in  fact  resolve 
isolated  shocks  as  well  as  an  exact  Riemann  solver.  Consequently,  using  this  robust 
solver  in  the  vicinity  of  strong  shock  waves  does  not  necessarily  pollute  a  scheme’s 
resolution,  as  would  be  the  case  if  artificial  dissipation  were  used  to  augment  the 
dissipation  provided  via  upwinding.  For  many  inviscid  calculations  the  amount  of 
pollution  proves  to  be  negligible,  and  whilst  some  degradation  would  be  expected  for 
the  case  of  a  strong  shock  interacting  with  a  boundary  layer,  it  may  well  be  unneces¬ 
sary  to  employ  the  HLLE  solver  in  such  a  situation  because  of  the  extra  dissipation 
provided  by  the  real  viscous  terms.  Secondly,  although  the  HLLE  switching  function 
requires  a  tunable  parameter  a,  the  retuning  of  this  parameter  is  less  involved  than 
the  retuning  of  an  artificial  dissipation  mechanism;  in  general,  it  is  far  simpler  to 
determine  where  extra  dissipation  should  be  added,  than  it  is  to  determine  how  much 
extra  dissipation  to  add.  For  many  problems,  assuming  that  shocks  are  located  as 
described  above,  a  sensible  threshold  on  the  shock  strength  can  be  specified  a  pri¬ 
ori.  Lastly,  it  should  be  noted  that  our  strategy  of  switching  Riemann  solvers  may 
not  prove  suitable  for  those  implicit  schemes  which  require  that  the  numerical  flux 
function  be  differentiable. 

Figure  7  shows  how  the  HLLE  solver  may  be  used  to  correct  the  tendency  of 
Roe’s  scheme  to  produce  kinked  Mach  stems,  c.f.  Figure  4.  For  this  calculation 
the  HLLE  switching  function  was  tuned  such  that  it  would  only  be  activated  by  the 
incident  shock,  and  the  principal  Mach  stem;  a  was  simply  set  to  half  the  strength  of 
the  incident  shock.  Note  that  apart  from  the  region  near  the  Mach  stem,  these  new 
results  are  very  similar  to  the  old  ones.  This  shows  that  the  HLLE  scheme  has  had  no 
adverse  affect  on  the  the  resolution  of  Roe’s  scheme.  Similarly,  Figure  9  shows  how  the 
carbuncle  phenomena  may  be  circumvented,  c.f.  Figure  3.  Here  we  have  restricted  the 
HLLE  solver  to  cells  near  the  stagnation  line  in  order  to  demonstrate  how  localized  the 
failing  of  Roe’s  scheme  really  is.  In  practice,  however,  we  would  advocate  using  the 
HLLE  scheme  along  the  whole  length  of  the  bow  shock,  so  as  to  maximize  robustness 
without  compromising  resolution.  Again,  a  sensible  value  of  a  can  be  found  a  priori 
by  using  some  large  fraction  of  the  shock  strength  along  the  stagnation  line  which 
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can  be  estimated,  given  the  freestream  Mach  number,  by  assuming  the  flow  is  locally 
one-dimensional.  As  shown  in  Figure  8,  the  HLLE  solver  may  also  be  used  to  good 
effect  to  prevent  Godunov’s  scheme  from  admitting  expansion  shocks,  c.f.  Figure  1. 
Here  we  have  employed  the  HLLE  solver  along  the  sonic  line,  and  in.  regions  where 
the  expansion  waves  are  strong. 

Having  presented  the  gist  of  our  strategy,  we  see  little  point  in  trying  to  sell 
a  particular  combination  of  solvers.  Starting  with  some  high  resolution  Riemann 
solver,  whose  choice  will  inevitably  be  a  matter  of  personal  taste,  the  correct  com¬ 
bination  of  solvers  will  depend  both  on  that  schemes  weaknesses  and  on  the  specific 
application  in  hand.  In  turn,  the  combination  of  Riemann  solvers  will  dictate  the 
choice  of  switching  functions.  Therefore,  we  shall  resist  the  temptation  to  recom¬ 
mend  a  specific  course  of  action,  instead,  we  present  two  simulations  that  show  how 
an  adaptive  Riemann  solver  might  be  used  to  good  effect.  Briefly,  both  simulations 
were  done  using  the  two-dimensional  analogue  of  the  one-dimensional  Euler  equations 
given  in  Section  3.  These  equations  were  integrated  using  the  two  step  finite-volume 
scheme  which  is  attributable  to  Hancock[25].  This  scheme  employs  van  Leer’s  MUSCL 
approach[26]  to  achieve  a  second-order  extension  to  Godunov’s  method,  hence  differ¬ 
ent  Riemann  solvers  may  be  slotted  directly  into  the  scheme  so  as  to  vary  the  flavour 
of  the  upwinding.  Although  the  calculations  were  performed  using  an  adaptive  mesh 
algorithm[17,  18],  the  mesh  refinement  monitor  function  was  such  that  the  calcula¬ 
tions  employed  a  nominally  uniform  cartesian  mesh. 

Our  first  example  concerns  the  simulation  of  a  strong  shock  wave  diffracting 
around  a  90°  corner,  the  shock  Mach  number  and  the  ratio  of  specific  heats  are 
5.09  and  1.4  respectively.  We  have  computed  this  test  problem  using  a  combination 
of  three  different  Riemann  solvers;  Toro’s  linearized  Riemann  solver  was  used  to  per¬ 
form  the  MUSCL  reconstruction  step  of  Hancock’s  scheme  as  described  by  Quirk[18], 
and  the  upwinding  step  was  performed  by  adaptively  selecting  between  Roe’s  solver 
and  the  HLLE  solver.  The  parameter  a  used  by  the  switching  function  (6)  was  set  to 
1  so  a.s  to  limit  the  HLLE  solver  to  the  incident  and  diffracted  shock  fronts,  and  to  a 
small  region  near  the  apex  of  the  corner.  Figure  10  shows  a  Schlieren-type  snapshot 
taken  from  this  simulation,  the  different  shades  of  grey  depict  the  magnitude  of  the 
gradient  of  the  density  field,  the  darker  the  shade  the  larger  the  magnitude;  details  of 
this  shading  procedure  are  given  in  Appendix  A.  Here,  it  is  not  our  intention  to  assess 
the  accuracy  of  these  results,  the  interested  reader  may  do  this  using  the  experimental 
results  of  Bazhenova  et  a/.[l],  and  the  computational  results  of  Hillier[7].  Instead  v'e 
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wish  to  illustrate  the  fact  that  certain  Riemann-solver  failings,  if  left  unaddressed, 
can  place  an  upper  limit  on  the  resolution  of  simulations  that  may  be  performed. 

Consider  the  consequences  of  disabling  the  HLLE  switching  function,  so  that  Roe’s 
solver  alone  is  used  for  Mie  upwinding  stage  of  Hancock’s  scheme.  The  tendency  of 
Roe’s  solver  to  allow  odd-even  decoupling  to  occur  along  a  planar  shock  wave  which 
is  aligned  with  the  grid  will,  sooner  or  later,  cause  this  simulation  to  come  to  grief, 
see  Figure  1 1,  thus  precluding  the  possibility  of  performing  very  detailed  simulations. 
By  way  of  comparing  the  resolution  of  these  two  sets  of  results,  for  Figure  10  there 
are  560  mesh  cells  from  the  apex  of  the  corner  to  the  point  where  the  Mach  stem 
meets  the  wall,  for  Figure  11  there  are  only  120  cells.  The  question  of  whether  or 
not  our  adaptive  mesh  algorithm  contains  some  flaw  which  exasperates  the  odd-even 
decoupling  is  largely  academic.  The  fact  remains,  that  Roe’s  solver  is  susceptible  to 
this  mode  of  failure  whereas  the  HLLF  solver  is  not.  Whether  the  initial  stimulus 
comes  from  a  distorted  mesh  as  in  Section  2  or  from  some  component  of  the  mesh 
adaption  scheme,  as  seems  likely  here,  is  immaterial. 

So  ''=■  no‘  to  leav^'  the  impression  that  the  above  shortcoming  is  somehow  pecu¬ 
liar  to  Roe’s  method,  we  present  a  second  set  of  results  which  are  taken  from  the 
interaction  of  a  planar  shock  wave  with  a  half-diamond;  the  shock  Mach  number  is 
2.85,  the  ratio  of  specific  heats  is  1.4,  and  the  angle  at  the  apex  of  the  diamond  is 
90°.  As  before,  we  have  run  this  test  problem  using  a  combination  of  three  different 
Riemann  solvers,  but  this  time  we  have  substituted  an  exact  Riemann  solver[24]  in 
place  of  Roe’s  linearized  solver.  Figure  12  shows  a  Schlieren-type  snapshot  from  this 
calculation,  note  that  some  800  cells  cover  the  width  of  the  diamond  so  this  calcula¬ 
tion  is  well  resolved.  Also,  as  an  aside  we  note  that  the  quality  of  these  results  may¬ 
be  gauged  by  comparing  them  with  the  experimental  results  given  by  Glass  ct  a/. [6]. 
Once  again,  if  the  HLLF  switching  function  is  disabled,  the  simulation  is  ruined  by 
the  odd-even  decoupling  that  occurs  along  the  length  of  the  incident  shock,  see  Fig¬ 
ure  13.  Note  that  this  second  calculation  is  of  lower  resolution  than  the  first,  only 
400  cells  cover  the  width  of  the  diamond. 

In  this  section  we  have  attempted  to  show  that  the  robustness  of  Godunov-type 
schemes  may  be  improved  by  employing  an  adaptive  Riemann  solver,  where  the  spe¬ 
cific  flavour  of  upwinding  is  altered  to  suit  the  local  flow  conditions.  If  used  sensibly, 
this  strategy  can  overcome  most  known  failings  of  individual  solvers.  Despite  our 
efforts,  we  recognize  that  the  majority  of  shock-capturing  practitioners  will  continiK' 
to  use  artificial  dissipation  as  a  band  aid  to  fix  a  particular  Riemann  solver  at  tlu' 
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first  signs  of  any  failing,  simply  because  it  is  expedient  to  do  so.  Whilst  vve  find  this 
approach  disappointine;,  the  principal  aim  of  this  paper  is  to  emphasiz<'  the  fact  that 
most  of  the  Riemann  solvers  that  are  in  common  use  must  be  augmented  in  some 
way,  if  they  are  to  be  used  for  the  purpose  of  producing,  genuinely,  high  resolution 
simulations  of  shock  hydrodynamic  phenomena. 


Figure  7:  The  HLLE  scheme  can  be  used  to  circumvent  the  tendency  of  Roe’s  method 
to  produce  kinked  Mach  stems:  (a)  Pressure  Contours,  (b)  HLLE  switching  function. 


(a)  (b) 


Figure  8:  The  HLLE  scheme  can  be  used  to  prevent  Godunov’s  method  from  produc¬ 
ing  expansion  shocks  (a)  Density  Contours,  (b)  HLLE  switching  bmction. 


Figure  9:  The  HLLE  scheme  can  be  used  to  circumvent  the  carbuncle  peh 
(a)  Density  Contours,  (b)  HLLE  switching  function. 


-  Xi  - 


Figure  10;  A  Scliliereii-type  snapshot  from  the  diffrartion  of  a  strong  shock  wave 
around  a  90“’  degree  corner. 


Figure  11:  On  its  own,  Roe’s  approximate  Riemann  solver  cannot  be  used  to  repro¬ 
duce  the  resolution  of  the  simulation  shown  in  Figure  10. 


-  2.5 


Figure  12:  A  Schlieren-type  snapshot  from  the  interaction  of  a  planar  shock  wave 
w'ith  a  half-diamond. 
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5  Conclusions 

Unless  the  dissipation  provided  via  upwinding  is  augmented  by  some  other  mecha¬ 
nism,  any  Godunov-type  scheme  built  upon  a  single  Riemann  solver  will  probably  be 
flawed.  For  example,  when  subjected  to  some  small  but  systematic  form  of  pertur¬ 
bation,  most  popular  Riemann  solvers  for  the  Euler  equations,  including  the  exact 
solver,  cannot  prevent  odd-even  decoupling  occurring  along  the  length  of  a  strong 
shock  wave  which  is  aligned  with  the  computational  mesh.  Thus  far,  this  flaw  has 
gone  largely  unnoticed  simply  because  it  is  only  exposed  by  very  high  resolution  sim¬ 
ulations.  However,  given  that  the  required  perturbations  can  arise  quite  innocuously, 
this  mode  of  failure  should  prove  fairly  widespread  once  genuinely  high  resolution 
simulations  become  commonplace  due  to  increases  in  computer  power. 

Although  most  flaws  can  be  controlled  by  the  judicious  use  of  a  small  amount 
of  artificial  dissipation,  to  do  so  necessarily  leads  to  a  reduction  in  the  resolution  of 
the  scheme.  We  favour  an  alternative  approach,  whereby  the  the  failings  of  any  one 
Riemann  solver  are  circumvented  by  combining  it  with  one  or  more  complementary 
solvers.  In  essence,  we  advocate  selecting  the  precise  flavour  of  upwinding  to  suit 
the  flow  data.  Admittedly,  this  synergetic  strategy  is  not  as  aesthetically  pleasing  as 
having  a  single  Riemann  solver  for  all  occasions,  but  we  have  shown  that  it  can  be 
made  to  work  quite  effectively.  Besides  which,  Riemann  solvers  are  sometimes  touted 
as  being  a  solution-adaptive  technique,  so  the  concept  of  an  adaptive  Riemann  solver 
is  not  that  contrived. 

Looking  to  the  future,  it  is  to  be  hoped  that  genuinely  multi-dimensional  Rie¬ 
mann  solvers  will  overcome  many  of  the  shortcomings  of  today’s  dimensionally-split 
schemes.  However,  given  the  way  in  which  the  present  shortcomings  have  been  stum¬ 
bled  across,  these  multi-dimensional  schemes  may  themselves  arrive  complete  with 
subtle  failings  with  which  to  ensnare  the  unwary. 
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A  Schlieren-type  Plots 


The  plots  shown  in  Figures  10  and  12  depict  the  magnitude  of  the  gradient  of  the 
density  field,  viz 


2 


and  hence  they  may  be  viewed  as  idealized  Schlieren  images.  So  as  to  accentuate 
weak  flow  features  the  following  nonlinear  shading  function  has  been  used, 


shade  =  exp{—kip) 


Here,  k  is  a  constant  and  4’  is  a  weighting  function  given  by 


|Vp|  -  |Vp|o 

|Vpli-lV;.|o’ 


where 


l^p|o  —  ^o|^pl7noa:i 


ko  and  k-[  being  constants.  Note  that  shade  is  limited  to  values  between  0  and  1, 
so  for  a  24  bit  colour  graphics  system  the  grey  level  shade  may  be  converted  to  an 
<  R,G,  B  >  triplet  using 


<  255  *  shade,  2b5  *  shade,  255  *  shade  >  . 

For  both  Figures  the  constants  k,  ko,  and  were  set  to  5,  0.05  and  —0.001  respec¬ 
tively. 


REPORT  DOCUMENTATION  PAGE 


Form  Approvor^ 
0MB  No  0 '^04-01 88 


1 


O-iO'i'  ’POC'T  ’'c  Di,''ce^  -  .p,  j., 1  o»’-  '“SO ;  t-*"  *  -  •  & .  -p.v  ^  :  r^st ' -  .•- 1 '  “  •  s'  -• :  a  ‘r.'  sox'-- 

gatr'ef"' J  iro  t''e  data  ne^aea.  ara  .  :~-o-i»!.r' :  -j''o  ■oa  i-*  -r-*- t^’-z  ;  ®s’  "•  .  ■  '  .  :-‘'r  t  '  s 

cOi'tf-rti'On  :  ‘  'c'  t  j  'ois  o-'3'="'  t  .'•  "Paao- .  .ces  '.  a'*  a*  c  Od#*' ■•'■■,'■-  J'‘o  d-  ’•'s  .  :  .e'*e'S- 

Davii  Miqhrtav.  So.te  ’20*1  **nir  qicn.  .  -i  22202 -4  302  a^'O  *.:  •’'<►0**-  ^  '*  ‘.‘  -i’  ja***^**''':  -''2  -'jC  }•?’  ^  ‘D<*r  ac'*  ^oc^i  t- r  :  'C-iO  '^^S)  /.  as^'  ''-:r  .:  2. :‘,Z: 

1.  AGENCY  USE  ONLY  (Le^ve  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AN 

November  i992  Contraf-Tor  R 

D  DATES  COVERED 

4.  TITLE  AND  SUBTITLE 

A  CONTRIBUTION  TO  THE  GREAT  RIEMANN  SOLVER  DEBATE 

S'.  FUNDING  NUMBERS 

C  NASl-18605 

C  NASl-19480 

WU  505-90-52-01 

6.  AUTHOR(S) 

James  J.  Quirk 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADORESS(ES) 

Institute  for  Computer  Applications  in  Science 
and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 

Hampton,  VA  23681-0001 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

ICASE  Report  No.  92-64 

9.  SPONSORING.  MONITORING  AGENCY  NAME(S)  AND  AODRESS(ES) 

National  Aeronautics  and  Space  Administration 

Langley  Research  Center 

Hampton,  VA  23681-0001 

10.  SPONSORING  /  MONITORING 

AGENCY  REPORT  NUMBER 

NASA  CR-191409 

ICASE  Report  No.  92-64 

11.  SUPPLEMENTARY  notes  Submitted  to  International 

Langley  Technical  Monitor:  Michael  F.  Card  Journal  for  Numerical 

Final  Report  Methods  in  Fluids 

124.  DISTRIBUTION ;  AVAILABILITY  STATEMENT 

Unclassified  -  Unlimited 

Subject  Category  02,  64 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (MaKimum  200  words) 

The  alms  of  this  paper  are  threefold:  to  increase  the  level  of  awareness  within  the 
shock  capturing  community  to  the  fact  that  many  Godunov-type  methods  contain  subtle 
flaws  that  can  cause  spurious  solutions  to  be  computed;  to  identify  one  mechanism 
that  might  thwart  attempts  to  produce  very  high  reolution  simulations;  and  to  proffer 
a  simple  strategy  for  overcoming  the  specific  failings  of  Individual  Rlemann  solvers. 

14.  SUBJECT  TERMS 

Rlemann  solvers;  shock  waves;  numerical  artifices 

15.  NUMBER  OF  PAGES 

33 

16  PRICE  CODE 

A03 

17.  SECURITY  CLASSIFICATION 

OF  REPORT 

Unclassified 

18  security  CLASSIFICATION 

OF  THIS  PAGE 

Unclassified 

19.  SECURITY  CLASSIFICATION 

OF  ABSTRACT 

20.  LIMITATION  OF  ABSTRACT 

MSN  7SA0  0’ -280-5500  Standard  form  298  ^Rev  2  89) 


P'r>^.  r,D»>g  Qy  A'A  ST3 
'.^2 

NASA  l.,thi;I.A 


